Monte Carlo with Absorbing Markov Chains: 
Fast Local Algorithms for Slow Dynamics 

M. A. Novotny 
Supercomputer Computations Research Institute, Florida State University, Tallahassee, FL 
J/-/. 32306-4052 

0\ ■ (February 6, 2008) 

i— s 

<3N ■ Abstract 

CO 

£> ■ A class of Monte Carlo algorithms which incorporate absorbing Markov chains 



OO 

o 



is presented. In a particular limit, the lowest-order of these algorithms re- 
duces to the n-fold way algorithm. These algorithms are applied to study 
ON ■ the escape from the metastable state in the two-dimensional square-lattice 



Ctf 



O 
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Monte Carlo (MC) methods |IJ have become indispensible tools for nonperturbative 
calculations in numerous fields, including materials science, high-energy physics, chemistry, 
biology, engineering, and economics. These methods are used for two fundamentally different 
purposes: to calculate time-independent quantities (statics) and to simulate time series 
(dynamics). In the former case, the slow relaxation observed, e.g. near phase transitions 
(critical slowing down) and at low temperatures, is merely a nuisance that has been overcome 
by a number of new MC algorithms, including cluster algorithms |§ , vertex algorithms |§ , 
multicanonical algorithms ||, and hybrid MC algorithms f|. Such algorithms can be many 
orders of magnitude faster than standard MC methods. However, they all replace the 
standard MC dynamic with a different dynamic. Consequently, although such algorithms 
may be very efficient in the calculation of static quantities, information about the kinetics 
of the original MC dynamic cannot be obtained. There are many instances where the 
kinetics, rather than just the statics, is of physical importance. Recently, MC methods using 
constrained cluster-flipping algorithms have been proposed in order to obtain information 
about the long- wavelength kinetics of a system |J. However, in such methods the local 
dynamic is modified, and universality arguments must be made to relate the results to the 
dynamic of the original system. 

In this Letter, a different class of accelerated MC algorithms is presented that does 
not change the original MC dynamic, and can therefore be used to simulate time series 
with very large separations in time scales. These algorithms incorporate into the standard 
MC algorithm absorbing Markov chains (AMC). The acronym MCAMC stands for Monte 
Carlo with absorbing Markov chains. With some additional approximations, the lowest- 
order MCAMC algorithm corresponds to the n-fold way algorithm 0. We demonstrate 
that MCAMC algorithms can be many orders of magnitude faster than either the standard 
MC or n-fold way algorithms. 

In a Markov process the probability of going from the current state of the system to 
another state depends only on these two states, not on the history of how the current state 
was reached. Standard MC algorithms are typical examples of Markov processes. A Markov 



process with a finite state space and a discrete time step is a Markov chain. The Markov 
chain is often written as a Markov matrix, M, the elements of which are the transition 
probabilities between the states. The time evolution of the probability distribution vector 
v 1 is then given by v^im + 1) = {/ r (r7i)M. An AMC is one in which at least one state has 
the property that transitions out of that state are forbidden. 

The Markov matrix associated with an AMC with q absorbing states and s transient 
states is given by || 

(i-qxq Oqxs \ 
(1) 
-T^s X q J- s X s / 

where I is the identity matrix, is the zero matrix, and the subscripts show the size of 
each submatrix. The elements of the transient submatrix, T, are the transition probabilities 
between the s transient states of M. The elements of the submatrix R are the transition 
probabilities from each of the s transient states to the q absorbing states. Starting from an 
initial vector vf, the probability of still being in any of the s transient states after m time 
steps is given by vfT m e, where e is the vector with all elements equal to 1. Although vf 
could be any normalized probability distribution over the transient states, in this work we 
will take vf to represent one particular transient state. It is easy to show that one can obtain 
the time m for exiting to one of the q absorbing states from the solution of the equation 

vfT m e<r<vfT m ~ 1 e, (2) 

where r is a random number uniformly distributed on the interval (0,1]. Note that the 
stochastic variable m is independent of which of the q states the AMC ends up in. The 
elements of the vector of q unnormalized probabilities 

Q T = tfTET.Xx, (3) 

give the probability of exiting to each corresponding state, given that the system has exited 
from the transient subspace in m time steps. After normalization, this probability distri- 
bution can be used to pick a particular state into which the AMC exits. This is done by 



using another uniformly-distributed random number. Eqs. (Q) and (0) are the only equations 
needed to include absorbing Markov chains within a Monte Carlo simulation. In practice, it 
is computationally sound policy to obtain all the eigenvalues and eigenvectors of T and uti- 
lize the spectral decomposition of T within the AMC portion of the algorithm to numerically 
solve Eqs. (g) and (§). 

The MCAMC algorithm has the following steps: 0) The system is in an initial state 
represented by iff. 1) Write an sxs transient matrix that includes at least the initial 
configuration, but ideally will include other states (which may represent many configurations 
of the system that are related by symmetry). The additional states included as transient 
states should be those the current configuration has the largest probability of exiting to. 
2) All configurations that the transient states can exit to in one time step must be included 
in the absorbing states. 3) Use Eq. (|2|) to find the time spent in the transient state subspace, 
and Eq. @ to find a new configuration for the system. 4). Iterate the procedure using the 
new configuration as the next initial configuration. Thus, once the system has exited to one 
of the q absorbing states, this decides the initial configuration of the system for a different 
AMC with different transient and absorbing subspaces, and the process is iterated. 

To illustrate the MCAMC algorithm, we apply it to the square-lattice nearest-neighbor 
Ising ferromagnet in a magnetic field. The Hamiltonian is given by 7i = —JJ2(i,j) s i s j ~ 
HJ2i s «) where the spins Sj=±l. The sums run over all nearest-neighbor pairs and over all 
N=L 2 sites, respectively. Periodic boundary conditions are used. The isotropic two-body 
coupling constant is given by J>0 and the magnetic field by H. To study the decay of a 
metastable phase, we apply a negative magnetic field at a temperature below the critical 
temperature T c , and start with the configuration of all spins +1 (which we call the C+ 
configuration). Standard droplet theory (for a recent review on metastability see Ref. [[|) 
shows that in order for the metastable phase C+ to decay, one or more critical droplets must 
nucleate via the random superposition of microscopic fluctuations. Since the critical droplet 
is of a certain size, the average time before its creation, and thus the metastable lifetime, 
is much longer than the microscopic timescales. Although any local dynamic can be used 



within the MCAMC framework, in this letter only Metropolis updates |l[ are performed. We 
measure the number of Monte Carlo steps per spin (MCSS) until a configuration is reached 
which has an equal number of +1 and —1 spins. The average lifetime, r, of the metastable 
state is found by averaging over a number of starts from the C+ configuration. 

A standard MC algorithm randomly chooses a spin, and then decides whether or not 
to flip it. Each spin in the LxL lattice, with periodic boundary conditions, is in one of 10 
possible energy classes which are determined by how the spin is oriented with respect to the 
applied field and the nearest-neighbor spins @. The number of spins in class i is n«, and 
the probability of flipping a spin in class i once it has been chosen is p(i). 

The s=l MCAMC algorithm is defined to have a single transient state, which is the cur- 
rent state of the spin configuration. When spins belonging to each class are grouped together 
the submatrix Ri x io = (nip(l), • • • , rii p(10))/N. Define Qo=0 and Qj=Y^i=i n iP(i)- Since 
M is a Markov matrix, its row-sums are unity, so the transient matrix is Ti x i=l— Q\q/N . 
The time increment to flip a spin is found from Eq. (Q) to be m>ln(r)/ln(l— Qio/N)>m— 1. 
If one relaxes the constraint that m is an integer, and expands to lowest order for small 
Qio/N, this equation becomes fn=— N\n(r)/Qi , and one has the standard n-fold way algo- 
rithm |7|]. In both the s=l MCAMC and the n-fold way algorithm, another random number 
is used to choose which of the 10 classes to pick via Qj-i<rQio<Qj, and a spin from class j 
is randomly picked and flipped. As shown in Fig. |], at low temperatures the s=l MCAMC 
algorithm can be many orders of magnitude faster than the standard MC algorithm. At 
J/T=3 the speed improves by a factor of about 10 7 . 

At low temperatures in the s=l MCAMC algorithm, once a spin has flipped from the 
C+ configuration, it is extremely probable that in the next iteration this spin will be chosen 
again to flip, and one returns to the C+ configuration. This problem can be addressed 
by going to the s=2 MCAMC algorithm. Whenever the starting configuration is the C+ 
configuration or a configuration with only one spin overturned from the C+ configuration, 
one includes these two states in the transient matrix T. The normal s=l MCAMC algorithm 
is used if the spin configuration is other than one of the transient states. As shown in Fig. 



|l|, for J/T=3 the speed improves by a factor of about 10 2 over the n-fold way algorithm. 

One can continue to increase the number of states in T. Figure [I] also includes results 
from s=3 MCAMC, where the three states in T were C+, C+ with one overturned spin, and 
C+ with two nearest-neighbor overturned spins. At J/T = 3 this s=3 MCAMC algorithm 
improves the speed by a factor of about 10 2 compared to the s=2 MCAMC algorithm. 

The MCAMC algorithm offers a further advantage, in that the statistical error in the 
average lifetime r is smaller than from standard MC. The average lifetime of an AMC 
that starts from state k is given by v£ Ne, where the fundamental matrix is defined by 
N=(I— T) _1 ||. For s>l, the initial configuration in which the AMC starts is C+ only 
once, since it must re-enter C+ from a configuration C+ with one overturned spin. Since 
the contribution to the lifetime that comes from the escape from C+ is exact, this reduces 
the total error in r. 

It has been shown that at low enough temperatures r is related to the height of the 
lowest energy barrier which must be reached by flipping one spin at a time starting from 
C+ [ID]. In the limit of low temperature and large L, the discreteness of the lattice gives 



an important contribution to r. Theorem 3 of Ref. [IDJ states that then 

k B T\n(L 2 r) = T(H, J) = 8J£ C - 2\H\(f c -4 + 1) (4) 

where £ C =\2J/\H\~\ and the notation \x] denotes the smallest integer larger than x. This 
result, which is restricted to 2J/\H\ not an integer and to \H\<4J, is shown in Fig. [I] 
as a solid line. It is in qualitative agreement with the measured values of r. A detailed 
comparison with predictions from Ref. [jlOj will be published elsewhere [11 . 



Using similar reasoning, it is possible to estimate at low temperatures the temperature 
dependence of the different MCAMC algorithms. By assuming that most of the CPU time 
is spent in the s—1 portion of the algorithm, at low enough temperatures the CPU time 
should be proportional to exp{[T(H, J)—T (H, J^/Zc^T}. Here To is 2J times the surface 
area minus 2\H\ times the volume of the largest compact lattice animal included in T. The 
dashed lines in Fig. |l] have these slopes and are drawn to go through the data point at the 
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lowest temperature available for a particular algorithm. 

Figure ^| shows values for r as a function of \H\ for two low temperatures. Note how well 
Eq. (^) fits the results. Eq. (|]) is only valid if the nucleation of a single droplet is responsible 
for flipping the lattice into the stable phase. A reasonable estimate for the crossover field 
(which has been called the dynamic spinodal field [O]) out of the single-droplet regime is 
if 1/2, the field at which the standard deviation of the lifetime is r/2 | 12|1 . Fig. |] shows that 



this is a reasonable estimate for the point at which the results deviate from Eq. (|J) at strong 
fields. 

At higher temperatures, the discreteness of the lattice becomes less important, and r 
should be given by the continuum droplet-theory prediction [|13|--|T6l 



ln(r) = E{T)/\H\ - (b + c) In \H\ + In (A(T)). (5) 

Here H(T) is related to the equilibrium zero-field surface tension, the equilibrium droplet 
shape [XT], and the spontaneous magnetization. A(T) is a non-universal function. Field- 
theoretical and numerical calculations give 6=1 |T^,|TB|,|TB[ . For dynamics described by a 
Fokker-Planck equation c=2 [|TJ|,|nj, and a recent MC study is consistent with b+c=3 fl2| |. 

Differentiating Eqs. (^) and (|5]) with respect to |-ff| _1 allows for direct comparisons with 
both the discrete droplet and continuum droplet predictions. This is shown in Fig. [3]. For 
T=0AJ a crossover is observed at weak \H\ to the continuum droplet prediction using b+c=3 
and the exact zero-field value for S(T) ||17|| . At weak \H\ Fig. |3| shows that the T=0.4J data 
are consistent with the theoretical predictions of Eq. (|5D with the theoretical values for S(T) 
and b+c. 

MCAMC algorithms can be utilized to study other interesting low-temperature effects 
in lattice-gas types of models. An example is the recent prediction that at low enough tem- 
peratures in the anisotropic Ising model nucleation occurs through square droplets, rather 
than through rectangular equilibrium Wulff droplets |19 . 



Additional applications of MCAMC algorithms may include studies of equilibrium prop- 
erties of systems where the n-fold way algorithm has been utilized. These include the 



anisotropic Ising model |20| (which is related via a Trotter-Suzuki decomposition to a quan- 
tum model in one lower dimension [pljj), and simulated annealing approaches to minimization 



However, the real strength of MCAMC algorithms lies in the study of slow dynamics 
in models with a limited number of local states. MCAMC algorithms will be particularly 
useful when there exists a small number of spins which are rapidly fluctuating, and major 
rearrangements of spins occur very infrequently. Examples would be the kinetics of spin- 



glass models [£f|, where the spins that fluctuate rapidly are the ones that find themselves 
in a local environment with zero energy difference between the possible spin orientations. 
Another example is phase-ordering kinetics in which a few particles undergo rapid random 
walks on long flat portions of interfaces, while the interfaces move extremely slowly. Similar 
reasoning also applies to simulations of molecular-beam epitaxial growth, where the n-fold 
way algorithm has been rediscovered at least twice |f23| , |24|| . In all such cases, using higher- 
order MCAMC algorithms should substantially decrease the CPU time required to obtain 
the static and dynamic information about the system without making any approximations 
to the dynamics. 
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FIGURES 
FIG. 1. The average CPU time for escape from the metastable state is shown as a function 

of the inverse temperature. This is for a 24x24 lattice in a field of |U|/J=0.75. For the standard 

Monte Carlo algorithm (x) the CPU time is directly proportional to the average lifetime r of 

the metastable state, and t (in units of Monte Carlo steps per spin) is plotted on the right-hand 

axis. All values of r were calculated by averaging over 10 3 escapes from the metastable state. The 

symbol x, with error estimates, for J/T>1 is calculated from the values of r from the MCAMC 

algorithms, while for J/T<1 the CPU time for a standard Monte Carlo algorithm is plotted. The 

CPU time from the ro-fold way algorithm, which corresponds to s=l MCAMC, is given by the 

symbol O. The timings for s=2 and s=3 MCAMC algorithms are given by the symbols + and 

□ , respectively. The vertical arrow marks the exact critical temperature. The solid line is the 

low-temperature discrete droplet result [pi for r with £ c =3 from Eq. (0). The dashed lines are 

estimates, described in the text, for the CPU times for the MCAMC algorithms. The horizontal 

arrows indicate that the solid line and symbols x are related to both vertical axes, whereas the 

other points and lines are only for the left-hand axis. The curvature near J/T~l is near the value 

where -ffi/2=3J/4 for this system size. Note the spectacular increases in speed with the MCAMC 

algorithms. 

FIG. 2. The average lifetime r, as a function of \H\, for escape from the metastable state is 
shown for a 24x24 lattice with T/J=0.2 (T/T c «0.088) (x) and T/J=0A (O)- Error estimates 
are plotted for each point. The solid lines are the low-temperature discrete droplet estimates for r 
from Eq. (||). The diagrams show the nucleating droplet, reading from the strongest fields with £ c 
of 1, 2, and 3. The dashed lines show the boundary between the various £ c regions. The vertical 
arrows mark the cross-over field Hi/ 2 [12], which separates the single-droplet regime (\H\<H]_i2) 



and the multidroplet regime (Hin<\H\) for the two temperatures for this lattice size. 
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FIG. 3. The temperature times the derivative of the logarithm of the average lifetime with 
respect to |-ff| _1 is shown as a function of \H\. This is for a 24x24 lattice with T=0.2J (x) and 
T=0.4J (O)- The solid curves are from the low-temperature predictions [1C] from Eq. (|j). Only 



results in the single droplet regime, \H\<Hi/ 2 , are shown. The dot-dashed and dashed horizontal 
lines correspond to the exact |l7j zero-field value for H(T)/T for T=0.2J and T=0.4J, respectively. 
The dashed inclined straight line includes the theoretical value b+c=3 in Eq. (|5|) for T=0.4J. 
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